Shear Viscosity of the Outer Crust of Neutron Stars: Ion Contribution 
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The shear viscosity of the crust might have a damping effect on the ampUtude of r-modes of 
rotating neutron stars. This damping has implications for the emission of gravitational waves. We 
calculate the contribution to the shear viscosity coming from the ions using both semi-analytical 
methods, that consider binary collisions, and Molecular Dynamics simulations. We compare these 
results with the contribution coming from electrons. We study how the shear viscosity depends on 
density for conditions of interest in neutron star envelopes and outer crusts. In the low density 
limit, we find good agreement between results of our molecular dynamics simulations and classical 
semi-analytic calculations. 
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I. INTRODUCTION 

Neutron stars are good resonators where oscillation 
modes can be excited. In particular, the r-modes of ro- 
tating neutron stars involve currents associated with very 
small density variations. These modes are unstable at all 
rates of rotation in a perfect fluid star [1]. The insta- 
bility is due to the emission of gravitational radiation, 
suggesting the possibility of gravitational wave detection 
with Laser Interferometer Gravitational- Wave Observa- 
tory (LIGO) [2J. This instability is expected to spin down 
newly born hot neutron stars [3]. However, the obser- 
vation of colder rapidly rotating neutron stars, suggests 
the existence of a damping mechanism of the r-mode in- 
stability. Several works have been done to explain this 
mechanism. For example, in Ref. |1] the damping mech- 
anism is suggested to be the result of a viscous layer at 
the interface between the solid crust and the fluid core. 
Other works discuss the r-mode dynamics of superfluid 
neutron stars [51 [^ finding that a core filled with neu- 
tron and proton superfluids limits the amplitude growth 
of the modes \J\. In addition, high multipolarity p-mode 
oscillations may impact the pulse shape of some radio 
pulsars f5]. For p-modes the primary restoring force is 
the pressure and the modes may be damped by the shear 
viscosity of the neutron star crust |^ . 

In the outer crust, the total shear viscosity has contri- 
butions from electrons and ions both of which transport 
momentum. Previous works have calculated the electron 
contribution to the shear viscosity in the Born approxi- 
mation |10j and including non-Born corrections [S]. Re- 
cently Horowitz and Berry calculated the electron con- 
tribution to the shear viscosity of non-spherical nuclear 
pasta phases [TT] . In this work we study the dependency 
of the shear viscosity with the density. We calculate the 
contribution of the ions to the shear viscosity in the neu- 
tron stars crust by two different methods: via Molecular 
Dynamics (MD) simulations and calculating momentum 
transport cross sections. The first method follows the 
Kubo formalism [12], and calculates the autocorrelation 
function of the pressure tensor. The second method con- 



siders binary collisions and allows one to consider both 
the classical and quantum systems. We focus our study 
to the case in which the ions form a dilute Hydrogen 
One Component Plasma (OCP). The conditions of tem- 
perature and density are chosen to reproduce those of 
the envelope and outer crust. However, the formalism 
can be applied to calculate other transport properties, 
such as diffusion coefficients. In particular the numerical 
procedure can be extended to study Multi-Component 
Plasmas (MCP). 

This paper is organized as follows: in Sec. |TT] we de- 
scribe our ion-ion interaction model, the Kubo formalism, 
and the semi-analytical procedure to calculate transport 
coefficients from transport cross sections. In Sec. |III| we 
present our numerical and semi-analytical results for the 
dependency of the shear viscosity with density. In Sec. 
IV we discuss quantum results vs classical ones. In Sec. 
V| we compare our results with the contribution coming 
from the electrons, and finally in Sec. |VI| we conclude. 
The appendix contains the description of the method em- 
ployed to calculate the phase shifts. 



II. FORMALISM 

A. Ion-Ion Interaction Model 

Electrons in the crust of a neutron star form a very 
degenerate relativistic gas that screens the interaction 
between ions. The ions form a plasma in the neutralizing 
electron background. It is convenient to characterize the 
strength of the ion interactions in terms of the Coulomb 
coupling parameter F, which is defined as the ratio of the 
average potential energy to the average kinetic energy, 



F = 






(1) 



where T is the temperature of the system, which we re- 
port in MeV (T[MeV] = fcsT), Z is the charge of the 



OCP, and 



Explicitly, the pressure tensor is given by 



"^ = (i) 



1/3 



(2) 



is the ion sphere radius. Finally, n is the ion density. The 
plasma frequency for a OCP is, 



M 



(3) 



where M is the mass of one ion. Quantum effects are 
expected to be important li T « Tp with Tp — LUp the 
plasma temperature. 

We describe the interaction between ions with a 
Yukawa potential. 



n»,j) 






-n,/A, 



(4) 



where 



is the distance between the zth and jth ions. 



Ae = TT^^^ /{2ekp) is the electron screening length with 
the electron Fermi momentum kp = (STr^rie)^'^, the elec- 
tron density is n^ = {Z)n, and e'^ ~ a the fine structure 
constant. In the OCP case all Zi correspond to the same 
ion species. 



B. Ion-ion Autocorrelation Function of the 
Pressure tensor 



We use the Kubo formalism to calculate transport co- 
efficients [12]. Kubo showed that linear transport coef- 
ficients L, could be calculated from a knowledge of the 
equilibrium fluctuations in the flux J associated with the 
particular transport coefficient. 



L = 



V 

fcRT 



dt{J{t)J{0)), 



(5) 



where V is the volume of the system, and T its tempera- 
ture. The integrand in Eq. (|5| is called the autocorrela- 
tion function. At zero time the autocorrelation function 
is the mean square value of the flux. At long times the 
flux J(t) at time t is uncorrelated with its value at t = 0, 
J(0) and the autocorrelation function decays to zero. 

In particular, to calculate the shear viscosity rj, we have 
made use of the autocorrelation function for the pressure 
tensor P, 



V 



SkeT 



I (l]^xy(i + io)Pxj;(io) Mi. (6) 



The average is taken over the three off-diagonal compo- 
nents {x,y), {y,z), (z,x), and over different initial times 
to- 



P. 



xyV-j — y 



^ rrijV^^ {t)vy^ {t) + -Y^ r^,, {t)Fy^^ (t) 



i^J 



(7) 

Here Vxi{t) is the x component of the velocity of the 
Jth ion at time t, r^-. is the distance in the x direction 
between the «th and jth ions, and Fy.. the y component 
of the force between them. 

This microscopic description allows us to calculate the 
viscosity of a gas, in which case we interpret the com- 
ponent Pij of the pressure tensor as the mean increase, 
per unit time and per unit area across a plane in the j 
direction, of the ith component of momentum of the gas 

m- 

Figure [T] shows the three different off-diagonal com- 
ponents of the pressure tensor. We need to calculate the 
correlation of this statistical noise. This is the aim of Eq. 
([6| . The components of the pressure tensor correspond 
to a simulation of OCP, for T = 1 MeV, n = 7.18 x 10"^ 
fm^^ , Z = 29.4, and A = 88. The total simulation time 
is 1.7 X 10^ fm/c and the time step is At = 50 fm/c. 




FIG. 1: (Color online.) The three off-diagonal components of 
the pressure tensor Eq. ([7|, corresponding to a simulation of 
OCP with Z = 29.4, A = 88, a temperature of 1 MeV and an 
ion density of 7.18 xlO"'' fm "^. 



In principle the integration in Eq. ^ will require an 
inflnite upper limit. In practice, we integrate up to a 
finite time tup- The choice of this upper limit requires 
caution. We want to find autocorrelations in the pres- 
sure tensor for long times, but at the same time if we let 
tup be too long the system becomes uncorrelated, and the 
statistical errors will increase. To solve this, we assume 
that after a perturbation the system relaxes exponen- 
tially. Then we find the relaxation time for the pressure 
tensor, and we integrate Eq. ([6| for up to few times this 
relaxation time. Figure [2] shows the result of integrating 
Eq. (Ig]), where the x axis represents different values of 
the upper limit t^p. The total simulation time is divided 
in four different runs and from them we obtain the error 



bars. Simulation parameters are as indicated in Fig. [T| 
The result for this case is 77 = 3.53(7) x 10"'^ fm "^ when 
tup — 5000 fm/c was chosen. 
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FIG. 2: (Color online.) Viscosity r; as a function of the upper 
time limit of Eq. mv. In this simulation all ions have Z = 
29.4, and A — 88. The temperature is 1 MeV and the ion 
density is 7.18 xlO"^ fm -^. 



C. Transport Cross Sections and Coefficients 



is a dimensionless momentum variable with a being the 
characteristic length scale of the potential. The quantity 
Ana^ with a — X^, for the Yukawa potential represents a 
normalizing cross section that renders the transport cross 
sections dimensionless. 

Since the interparticle distance limits the number of 
partial waves for scattering in a two-body system, we also 
introduce density dependent quantum transport cross 
sections by limiting the number of terms in the summa- 
tion to 



in 



[kn 



-1/3 



1/2] 



(11) 



where n is the number density and the quantity [c] de- 
notes the integer part of c. 

If the particles possess spin s, then the properly sym- 
metrized forms are: 



«( 



(") 






2s+l^Bose ^ 

s + 1 „(") 
2s+l "-Fermi 



_ (") 

^1 ^Fermi ' 



_ (n) 
-iIBos. 



for integer s. 



for half-integer s. 



(12) 

From Eqs. (l9| and (10 1, we note that phase shifts are 
the central physical input for the transport cross sections. 
Details of the calculation of the phase shifts, particularly 
for collisions of nuclei with large atomic number Z which 
involve a large number of partial waves and thus a large 
number of scattering phase shifts, are provided in the 
appendix. 



We obtain semi-analytical transport coefficients, for 
both quantum and classical systems, by calculating the 
appropriate cross sections. Our calculation follows the 
Chapman-Enskog formalism described in Ref. [13]. The 
quantity of interest is the transport cross section, 



(") 



2n 



dcos6'(l-cos"6i) 



da{k,e) 



dn 



(8) 



where the scattering angle and the coUisional differen- 
tial cross-section "^Ip,' I are calculated in the cen- 

"ii \c.m. 

ter of mass reference frame of the two colliding parti- 
cles with momentum Hk. For indistinguishable parti- 
cles, an expansion of the cross-section in partial waves 
X;i(2/-hl)(e*2*' -l)P,(cos6') and the orthogonality of the 
Legendre polynomials Pi simplifies the integrals above to 
the infinite sums 



^(1) = 



■,(2) = 



(2) 

47ra2 



,^'(2/+l)sin2(5,(x)). 



(9) 



2 , (^+l)(/ + 2) . 
^Z^ /o; I o\ ^™ [(>l+2[Xj - di{x)), 



(2/ + 3) 



(10) 



where the prime on the summation sign indicates the use 
of even I for Bosons and odd I for Fermions; x — ka 



The classical transport cross section is given by 



q(n) = 



1 

4a2 



(l-cos"(6l))d62, 



(13) 



with b the impact parameter, and 



e{b,E) 



200 



bdr 



'•o{b,E) ^2 



&2 



Vjr) 
E 



(14) 

the deflection angle which is a function of the energy 
E, b, and V{r). The quantity ro{b,E) is the distance 
of closest approach. We also study the case in which 
the upper limit for b in the integral above is set equal 
to the average interparticle distance « 7i~3. Note that 
this procedure introduces a density dependence into the 
classical transport cross sections. 

In the dilute limit, n ^ 0, the classical and quantum 
transport cross sections (7(2) fQj- systems with A = 88 
and A = 1 are shown in Fig. [3J As the energy tends to 
zero (or x —^ 0), the quantum cross section stays finite 
whereas the classical one diverges. For this reason, as we 
will see, the quantum result for shear viscosity will differ 
from the classical result only for very small temperatures. 
For increasing x, the summation over I in Eq. ( |10p in the 
quantum case leads to results that coincide with those of 



0.03 , 
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FIG. 3: (Color online.) Qua ntum and classical transport cross 
sections g'^' from Eqs. ( 10 1 and ( 13 1 at low energies for sys- 
tems with A — 1 and A — 88 in the dilute limit (n ^ 0). 



the classical case obtained using the deflection function 
(see Eq. ( 13 1). For the system with A = 88, keeping only 

5) is adequate in the energy 



five partial waves {1 = 
region where x < 0.25. 

The transport coefficients can be calculated from the 
omega-integrals 



i'"'*)(r)= / d7e-"'''72*+3g(")(x), 



(15) 



where 7 = , ?; , T is the temperature, and /i is the 
reduced mass. The formalism above includes only binary 
collisions. Therefore, we expect it would be accurate only 
at low densities. 

In the first order (of deviations from the equilibrium 
distribution function) approximation, the shear viscosity 
is given by [14j : 



Ml 



\Ht)Jj^^'\t) 



(16) 



where X{T) = h/^/2TrMkBT is the thermal de-Broglie 
wavelength. The quantity 



V 



16^2 a^ 



(17) 



is a characteristic viscosity of the Yukawa potential with 
a = Xe- For later use, we define here the characteristic 
temperature fc^f = f?g (or X = (|)2). 



Equation (16) shows clearly that if 



,(2,2) 



T- 



independent (as for rigid-spheres with a constant cross 
section), the shear viscosity exhibits a T^/^ dependence 
which arises solely from its inverse dependence with A(T). 
For energy-dependent cross sections, however, the tem- 
perature dependence of the viscosity is sensitive also to 
the temperature dependence of the omega-integral. 
In the second order approximation [15_ 



M2 



(l + J^(T))(l±nA%(r)), (18) 



where ± means plus sign for Bosc and minus for Fermi 
statistics, n is the number density. 



,(2,2) 



3(7cjr - 2lu\ 



(2,3)s2 



2 LP'^) (77u;p-^) -t- 6^^^)) - 6 (cP))' 



and 



2-7/2 



(2,2) 
128^4/3 
33/2 ,,(2,2) 



(19) 



(20) 



It is worthwhile to note that at the first order of de- 
viations from the equilibrium distribution function, the 
viscosity is independent of density, unless density depen- 
dent cut-offs are used to delimit the quantum or classical 
transport cross sections. An explicit density dependence 
arises only at the second order. 



III. MD SIMULATIONS OF A DILUTE PLASMA 

We calculate the shear viscosity for a dilute plasma 
following the procedure described in Sees. |II A| and II B 



The results presented here correspond to Hydrogen, Z 
1 and A = 1. The box length of the simulation volume 
for the parameters used is L ^ 100 fm. To minimize 
finite size effects in our MD simulations, it is needed that 
L 3> Ag , where Ag is the screening length of the Yukawa 
interaction, see Eq. (Bl. To assure this condition, we 
arbitrarily fix the electron screening length Ag to 10 fm. 
Also, we choose values of density and temperature such 
that the system is weakly coupled, F < 1. 

In tablclTlwe summarize the values of the ion density n, 
Coulomb parameter F, and simulation times used. The 
time step AtMD is 10 fm/c. We will use the results from 
these MD simulations here and later in Sec. llVl to com- 
pare with results obtained from calculations performed 
using the material in Sec. 



lie 



Here, we point out that our numerical simulations in 
the limit when n —^ are involved. To follow the trajec- 
tories of a highly dilute gas implies very small simulation 
time steps, and long simulation times in order to find 
correlations between the ions. 

Table |lT] shows results for the viscosity for different 
values of temperature and density. Again we use A = 1, 



TABLE I: Simulation runs at T = 0.1 MeV for a H dilute 
plasma. Nion is the total immber of ions in the system. The 
time step in the MD simulation is AtMD = 10 fm/c. Tw is 
the warm up time, Tm the measurement time, and Atp^^^ is 
the time step used to calculate the pressure tensor. 



^fm-^) Tvy(105fm/c) rM(10^fm/c) Atp^ 

8.0 10 

21.9 97.0 100 

0.2 100 100 

0.2 97.0 100 

47.2 12.5 10 

1.3 10 10 



4 T 



Nion 


r n(xlO 


500 


0.5 0.01 


100 


1.07 0.1 


100 


1.46 0.25 


100 


1.84 0.5 


500 


2.32 1.0 


500 


3.15 2.5 






O- 



£-2 
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quantum 





6 \ ..••! 



classical 



A = Z = 1 
a = 10 fm 



-2-10 1 

logio (kfiT, MeV) 



FIG. 
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4: (Color online. 



TABLE 11: (Color online.) Viscosity results (given by setting 
ft = 1) for a H dilute plasma and parameters used in the 
simulations. Nion, the total number of ions in the system, is 
500. The time step in the MD simulation is AImd, whereas 
Tmd is the total simulation time, and Atp^^^ is the time step 
used to calculate the pressure tensor. 



T n 


r r, 


Tmd 


^iP.y. 


AtMD 


(MeV) (fm-^) 


(fm-«) 


(lO'^fm/c) 




(fm/c) 


0.5 10"=^ 


0.46 3.9(3) X 10-^ 


7.35 


100 


5 


0.1 10"=^ 


2.32 3.28(7) x 10^'' 


1.72 


10 


10 


0.05 10-^ 


4.64 1.51(3) X 10"'' 


1.08 


10 


10 


0.01 lo-'' 


5.0 7.0(2) X lO-'' 


0.8 


10 


10 


0.001 10^5 


50 9.1(5) X 10^^ 


0.8 


10 


10 



and Z — \. The three first rows correspond to cases in 
which the density was kept constant. The viscosity is 
a monotonically increasing function of T. We find from 
our simulations that the dilute plasma behaves as a gas 
where higher temperatures provide higher momenta, and 
hence lead to larger momentum fiux. 

Table |ll] also shows the feature of increasing viscosi- 
ties for large densities. As the temperature decreases 
the differences due to changes in density are more evi- 
dent. We expect that as the temperature decreases the 
changes in the dependency of 77 vs n to be more dramatic 
making the plot steeper due to the effect of correlations 
between ions. For example, rj remains in the same order 
of magnitude ^ 10~^, when n changes by two orders of 
magnitude. On the other hand, comparing the runs at 
T = 0.05 MeV and T = 0.01 MeV (close values) in Table 
im we observe a change in rj of two orders of magnitude, 
when n changes in the same way. This is due to the fact 
that there is a larger change in the inter-ion distance, 
which increases the correlations between ions, than the 
change in the ratio of thermal to Coulomb energy F. 



Quantum and classical results from 
(16 1 for the first order shear viscosity [77]! for a H di- 
lute plasma. (Viscosities in this paper are given by setting 
h = 1). The lower-most solid curve is the classical dilute gas 
limit. The dots show the corresponding limit in the quantum 
case. The upward-bending solid (dashed) curves are results 
obtained using density dependent cut-offs in the quantum 
(classical) case and signal the onset of more than two-body 
effects. 



IV. DILUTE LIMIT QUANTUM AND 
CLASSICAL VISCOSITIES 



In this section we compare results based on the 
Chapman-Enskog formalism outlined in Sec 
those from the MD calculations in Sec. |III[ In Fig 
the quantum and classsical results from Eq, 
first order shear viscosity [ry]i for the system with A 
are shown as functions of temperature. As noted earlier, 
both classical (the lower-most solid curve) and quantum 
(the dotted curve) viscosities are independent of density 
at first order. As expected, the quantum results approach 
the classical ones at high temperature. With decreasing 
temperature, however, the viscosity for the quantum case 
is significantly larger than its classical counterpart due to 
the smaller transport cross sections (see Fig.[3|. 

The upward-bending solid (quantum case) and dashed 
(classical case) curves in this figure show results obtained 
by imposing the density dependent cut-offs on the angu- 
lar momentum in the quantum case, Eq. (Ill, and the 
impact parameter in the classical case. 



nqwith 

i 



(16|) for the 
1 



These results al- 
low us to establish the ranges of density and temperature 
for which the first order result is valid. For example, at 
a density of 10~^ fm~ , the quantum results for temper- 
atures below r '--^ 0.1 MeV are susceptible to more than 
two-body effects (not considered in this treatment) so 
that the first order result is not reliable. For the range of 
temperatures shown in this figure, effects of the density 
dependent cut-offs are more significant for the quantum 
case than for the classical case. 



ks T = . 1 MeV 
s = 1/2 




while to mention that when these corrections are signifi- 
cant, many-body correlations not considered here will be 
important. This physical effect is amply demonstrated 
by the results of the classical MD calculations which lie 
between the classical values of [77] i and [77] 2 • 



logio (n, fm" ) 



FIG. 5: (Color online.) Quantum and classical viscosities as 
a function of the ion density n, for a H plasma. The tempera- 
ture is 0.1 MeV. First and second order results from Eqs. ( |16[ ) 
and (18 1 are denoted by [ri]i and [ri]2, respectively. The filled 
squares are results of MD simulations. 



For the system with A — I, results of [q] up to second 
order from Eqs. (16 1 and (18 1 are shown in Fig. p^ as 
functions of density. The temperature in this case is 0.1 
MeV. The upward-bending curves show effects of density 
dependent cut-offs. Results from the MD calculations 
are shown as filled squares. Noteworthy features in this 
figure are: (1) the large differences between the quantum 
and classical results, and (2) the disposition of the MD 
results with respect to the classical results. These results 
also underscore the importance of incorporating quantum 
effects in MD calculations, even at low densities. 

In order to understand these results, we begin by not- 
ing that the characteristic shear viscosity and tempera- 
ture in this case are i) = 2.2 x 10^* fm^ (we have set 
h= 1) and fcsf = 2.6 MeV. For ksT = 0.1 MeV and 
A — 1, we find [r]]i (classical) — 1.72 x lO^'' fm~^ and 
[7j]i (quantum) = 1.82 x 10"^ fm"^. 

The fact that the first order viscosity in the quantum 
case is nearly ten times larger than the classical result can 
be understood by examining the ratio of the correspond- 
ing omega-integrals in Eq. ( 15 1. In the quantum case, the 
integrand of the omega-integral peaks at 7 ~ 2 for which 

we find Xpeak — 2\/2'kT/T ^1. At this value of Xpeak, 
the quantum transport cross section is nearly ten times 
smaller than the classical one (see Fig. Isl) which quanti- 
tatively accounts for the ratio of ~ 10 being sought. In 
physical terms, the differences in viscosities stem from 
the differences in the transport cross sections: in the 
quantum case, the cross section saturates at low energies 
whereas in the classical case the cross section diverges. 

The densities at which the viscosities begin to be de- 
pendent on angular momentum or impact parameter cut- 
offs are also evident from Fig. [s] In the classical case, 
effects of the density dependent cut-off enter at much 
higher densities than do differences stemming from going 
to a higher order (compare [ryji and [77] 2). It is worth- 




FIG. 6: (Color online.) Quantum and classical viscosities as 
a function of T/Tp for a H plasma. The plasma temperature 
Tp is given in Eq. ([3|. The densities are as indicated. First 
and second order results from Eqs. ( |16[ ) and ( |18[ ) are denoted 
by [r;]i and [77] 2, respectively. 



In order to appreciate how and when quantum ef- 
fects become important in a Hydrogen plasma, we show 
in Fig. [6] results of the dilute limit quantum and clas- 
sical first and second order viscosities for two densi- 
ties as functions of the ratio T/Tp, where Tp = Wp ~ 
27.5 (Z'^n/A)^/'^ MeV (with n in fm'^) is the plasma 
temperature in Eq. ([3]). Results obtained with density- 
dependent cut-offs are also included in this figure to illus- 
trate when three-particle effects become important. Re- 
sults for the intermediate densities (not shown in the fig- 
ure for the sake of clarity) are quantitatively different, 
but qualitatively similar. The filled circles are the re- 
sults of MD simulations corresponding to a temperature 
of 0.1 MeV at the densities shown in Fig. [5] The tra- 
jectory traced by the MD results is due to the density 
dependence of Tp . 

The upward shift in the magnitudes of the viscosities 
(which are intrinsically independent of density) with in- 
creasing density in Fig. ^ is caused by the n^' ^ depen- 
dence of the plasma temperature Tp. Quantum effects 
lead to viscosities that are significantly larger than the 
classical results as T/Tp decreases; the lower the densi- 
ties the larger are the differences. Put differently, the 
values of T/Tp for which the quantum and classical re- 
sults merge together increase as the density decreases. 
It must be borne in mind, however, that although the 
quantum effects are captured fully in the cross sections 
the thermal weightings are classical in our treatment 
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7; (Color online.) Same as Fig. |4] but for the system 
A = 88. 




here. With increasing density, effects of the Pauh prin- 
ciple and possible in-medium and many-body effects will 
become important. As shown by the second-order vis- 
cosities and the results obtained using density-dependent 
cut-offs, more than two-body effects become important 
when logio(T/Tp) ~ (2) for n = lO^^ (IQ-^) fm^^ 

The result for shear viscosity at first order is shown 
in Fig. |7]for the system with A=88. As the density de- 
creases, the classical (lines) and quantum (dots) results 
merge together as expected in the dilute limit. Differ- 
ences between the two cases occur only as the temper- 
ature decreases to values below the temperature corre- 
sponding to the plasma temperature. The character- 
istic shear viscosity and temperature in this case are 
r) ^ 1.2 X 10^5 fm"^ and fcsf = 4.4 x 10"^ MeV. For 
ksT = 1 MeV and A = 88, we find [r]]i (classical) = 
4.63x10-5 fm"^ and [r?]i (quantum) = 4.62 xlQ-^ fm'^ 

Results of [r]] up to second order are shown in Fig. [s] 
for the system with A = 88 as functions of density at a 
temperature of 1 MeV. There is virtually no difference 
between the classical and quantum results because the 
system is essentially classical. Effects of the density de- 
pendent cut-off enter at much lower densities than do 
differences stemming from going to a higher order (com- 
pare [77] 1 and [77] 2). 



Figure |9] shows how the viscosities of the heavy-ion 
plasma with A — 88 depend on the ratio T/Tp. Ex- 
cept for T/Tp <C 1, the classical and quantum results are 
indistinguishable. The values of T/Tp for which more 
than two-body effects are important are easily discerned 
from this graph by inspecting the results obtained using 
density-dependent cut-offs. 
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FIG. 8: (Color online.) Same as Fig. [sl but for the system 
with A = 88. 
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plasma with A=88. 



Comparison of Dilute Limit Viscosities with 
MD Results 



We find that MD results for the Hydrogen plasma are 
in good agreement with the classical semi- analytical cal- 
culations at low densities. On the other hand, the first- 
order semi-analytical solution predicts a constant behav- 
ior of the viscosity for all values of the ion density. MD 
simulations have taken into account correlations between 
ions which have not been included in the semi-analytical 
results. These correlations have a stronger effect as den- 
sity increases leading to a larger difference from the semi- 
analytical description. From our results we conclude that 
the viscosity for a dilute H plasma is constant at low den- 
sities and rapidly increases for higher densities. Results 
for T <C Tp, with Tp the plasma frequency, could involve 
quantum corrections and we expect our MD results to be 
inaccurate in that regime. 

We emphasize the agreement between our MD simu- 
lations and semi-analytic results in the low density clas- 
sical limit. This provides a significant check for both 



approaches. As the density increases, MD simulations 
fully include correlations between ions and are there- 
fore directly applicable at any density. In contrast the 
semi-analytic approach only works at low densities. The 
Yukawa interaction is reasonably long ranged. It can 
still be important at distances of a few or more screen- 
ing lengths Ae. Furthermore, there are a large number 
of other ions to interact with at large distances. This 
limits the applicability of the semi-analytic approach to 
low densities. 

On the other hand, the quantum semi- analytical vis- 
cosities for the Hydrogen plasma differ from those of clas- 
sical semi-anlytical and MD simulations by about an or- 
der of magnitude. This large difference points to the need 
for including quantum effects in MD simulations of light 
ions. However, see results in Sec. IV B in which density 



dependent screening lengths are used. 

The fact that the MD result rj = 3.53 x 10"^ fm "^ 
for the system with A = 88 at ksT = 1 MeV and 
71 = 7.18 X 10^^ fm^'^ is about two orders of magnitude 
larger than the semi-analyical result of 4.63 x 10~^ fm~ 
deserves some comments. From Fig. |7] we note that at 
this density the cut-off dependence (which signals the 
influence of more than two-body effects) in the semi- 
analtical viscosities sets in at a temperature of about 10 
MeV. Alternatively, Fig. |8] shows that at a temperature 
of 1 MeV, the dilute limit is reached only below a den- 
sity of 10^^ fm^ . Thus for the MD results to approach 
the semi-analytical results, either an order of magnitude 
higher temperature is required at this density or an order 
of magnitude lower density is needed at this temperature. 
For the values of density and temperature chosen, many- 
body effects will be important. Note that the second- 
order semi-analytical approach is an attempt, but only 
in a perturbative fashion (as an expansion in the param- 
eter (na^)). For this perturbative result to be trustable, 
it should not exceed the first-order result substantially. 
The MD calculations take into account many-body ef- 
fects to all orders, albeit classically in our calculations. 
For heavy-ion systems in the dilute limit, quantum ef- 
fects are not expected to play an important role because 
a lot of partial waves contribute. 
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FIG. 10: (Color online.) Viscosities for density dependent 
screening lengths as functions of density. 



the density dependent screening length Xe{n). As the 
screening length is density dependent, the shear viscos- 
ity exhibits distinct features as a function of density. In 
Fig. [To] results of viscosities for the Hydrogen plasma are 
shown as a function of density. Even at first order the 
shear viscosity is density dependent through the poten- 
tial and increases with increasing n. Corrections arising 
from second order deviations from the equilibrium distri- 
bution function contribute significantly even at low den- 
sities. Results obtained by imposing density dependent 
cut-offs indicate that more than two-body physics plays 
an important role at low densities in both the classical 
and quantum cases. The large screening lengths at low 
densities induce a large number (up to 1000) of partial 
waves to contribute in the quantum case. For moderate 
densities, however, significantly lower number of partial 
waves are needed to obtain convergent results. Both the 
classical and quantum results are nearly the same at low 
densities. However, the quantum results are larger than 
the classical ones as the density increases because the 
classical cross section is larger for a smaller screening 
length at a fixed energy (x — ka gets smaller). Note that 
imposing density cut-offs on the classical results mimics 
the quantum results without cut-offs. 



B. Dilute Limit Viscosities for Density Dependent 
Screening Lengths 

The electron screening in length in a charge neutral 
plasma is given by 



A. 



2efcp 



3.35 



1 



{Z)n 



1/3 



(21) 



For example, at a density of 10^^ fm^'', Ae — 336 (108) 
fm for (Z) = 1 (29.4). As MD calculations are time- 
wise prohibitive for very large screening lengths, we re- 
port here on the dilute limit semi-analytic calculations 
of Sec. |IIC| in which the parameter a in the Yukawa po- 
tential describing the ion-ion interaction is set equal to 



Additional insight can be gained by examining viscosi- 
ties as functions of T/Tp (see Fig. 111. Comparing the 
results for the two densities shown in this figure, we learn 
that the quantum and classical results begin to differ 
from each other at progressively lower values of T/Tp as 
the density decreases. Whereas the first order classical 
and quantum results differ by small amounts, the differ- 
ences grow with increasing density. Many body effects 
gauged through cut-offs appear together with second or- 
der corrections. 

An important feature that emerges from the above re- 
sults is that the extent to which the quantum results 
differ from the classical results is not as large as in the 
case when a was fixed for all densities as in Sees. IIIII and 
IV The feedback offered by the density dependence of 
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effects of cut-offs. Note that the effects of statistics are 
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FIG. 11: (Color online.) Viscosities for density dependent 
screening lengths as functions of T/Tp for the Hydrogen 
plasma (spin — 1/2). 
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FIG. 12: (Color online.) Same as Fig. |11[ but for the system 
with A = 88 (spin = 0). 



the screening length (this is in fact a many-body effect) 
serves to reduce the differences between the quantum and 
classical viscosities. 

In general, quantum effects are important for light ions 
and for small screening lengths Ae. In the limit Ae — + oo 
the potential reduces to a 1/r Coulomb potential for 
which the classical and quantum differential cross sec- 
tions agree. For Hydrogen with Z = 1, the screening 
length from Eq. ( pTj ) can be much larger than the fixed 
Ae = 10 fm used in the previous sections. This larger 
screening length leads to smaller quantum effects. (Note 
that the MD simulations are easier for Ag = 10 fm.) We 
conclude that quantum effects are probably small for re- 
alistic screening lengths. 



Figure [12] shows results of viscosities as functions of 
T/Tp for the system with A = 88 at two different den- 
sities. Whereas the first order classical and quantum re- 
sults agree over a wide range of T/Tp, the role of many- 
body effects studied through the influence of cut-offs is 
evident even for T/Tp ^ 1 . Effects of second-order cor- 
rections to viscosity enter at much lower T/Tp's than the 



opposite for spin 1/2 (Fig. Ill and spin systems 



V. COMPARISON WITH THE ELECTRON 
CONTRIBUTION 

The shear viscosity of the outer crust of a neutron star 
is determined by the contribution of the various matter 
components. In this way the total shear viscosity rjtot is 
given by the sum of the viscosity due to electrons r/g and 
the ions rj, rjtot — Ve + V- I^i this section we calculate the 
electron contribution for a single case, and compare with 
the ion contribution obtained in the previous section. We 
follow the formalism used by Chugunov and Yakovlev 
[S] . The shear viscosity coming from the electrons can be 
calculated from 



Ve 



riekpVF 



(22) 



where z/g is the effective electron collision frequency and 
vp (« 1) is the electron Fermi velocity. For dense matter 
this collision frequency is given by 

J^e = Vei + Uimp + l^ee, (23) 

where z/gi, J^imp, and Vee correspond to electron scattering 
from ions, impurities, and electrons, respectively. We 
assume electron-ion scattering is the dominant process. 
The electron-ion collision frequency is given by 
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where K^i is the effective Coulomb logarithm. 
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(24) 



Sn{q)dq, 



(25) 
where q is the momentum trasnfer. The lower limit go 
is in a liquid phase and go = {Qii^n)^'^ in a crystal 
phase [IB], the effective electron mass is m* — ep — 
{kp + ml)^''^ with kp the electron Fermi momentum and 
rrie the electron mass. Sjj^q) is the structure factor de- 
scribing electron-ion scattering, u{q) is the Coulomb in- 
teraction between an electron and a nucleus, and e{q, 0) 
is the dielectric function due to degenerate relativistic 
electrons, [TTIITH]. 

We calculate S',,(g) directly as a density-density corre- 
lation function using trajectories from our MD simula- 
tions, 



^.(q) - (p*(q)p(q)> 

Here the ion density p(q) is. 



l(p(q)>l' 



p(q) = 
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FIG. 13: (Color online.) Structure factor from MD simula- 
tions. The temperature is 0.1 MeV and n = 2.5 x 10~^ fm~^. 
The dotted curve is an anlytical fit from Ref. [16| . 



Ref. [9j calculates the electron contribution to the 
shear viscosity using the analytical fit for the effective 
structure factor proposed in Ref. [IS] plus some correc- 
tions included to take into account the form factor of 
the nucleus. As shown in [T2] our MD results for the 
Coulomb logarithms are comparable to the analytical re- 
sults of Ref. jl6j . The structure factor used to find the 
Coulomb logarithms is shown in Fig. [13] along with the 
anlytical fit from Ref. [TH]. The lowest value of q at 
which we can calculate Sri{q) is limited by go = 27r/L. 
We approximate S',,(q) by S',,(qo) for the low q region 
from to qQ. This value agrees with the Debye-Hiickel 
approximation at q = 0. 



Table [Hi] shows the viscosity obtained for a simulation 
with 500 ions at a temperature of 0.1 MeV and an ion 
density of 2.5 x 10~^ fm~^. Even with large quantum 
corrections, the contribution from the ion-ion scattering 
is small by orders of magnitude compared with the con- 
tribution from electron-ion scattering. Therefore, for this 
regime of densities and temperatures the ion contribution 
to the shear viscosity of the outer crust in a neutron star 
is negligible. Nevertheless, the use of MD simulations has 
allowed us to calculate both contributions. 

Our MD formalism can also be used to calculate the 
thermal conductivity n |19j . Again we expect the elec- 
tron contribution to k to dominate over the ion contri- 
bution at low magnetic fields. However, large magnetic 
fields suppress the electron contribution in directions per- 
pendicular to the field; then the ion contribution can be 
important. For example, Potekhin and Yakovlev [20' find 
that magnetic fields larger than 10^^ Gauss can lead to 
anisotropic temperature distributions. Therefore we ex- 
pect a field of 10^^ Gauss to start modifying electron con- 
tributions to the shear viscosity. However, it may take 
significantly stronger fields before the ion contribution to 
the viscosity dominates that from the electrons. 



TABLE III: Contributions to the shear viscosity rj (for the Hy- 
drogen plasma), coming from electron-ion scattering rje, and 
from ion-ion correlations rj. The system is at a temperature 
of 0.1 MeV and the ion density is 2.5 x 10"^ fm"^ 



Ae, (fm ^ 
0.27 



77e(fm ^) 7?(fm" 



26.61 



4.81 X 10" 



VI. SUMMARY AND CONCLUSIONS 

We have calculated the ion contribution to the shear 
viscosity in the outer crust of neutron stars. We use two 
different methods for this purpose. One of these is based 
on MD simulations while the other is semi- analytical and 
calculates the transport cross sections for classical and 
quantum systems. We find good agreement between the 
two methods in the low density, classical limit. 

Using the MD trajectories, we have used autocorrela- 
tion functions of the pressure tensor in a Hydrogen OCP 
to calculate its viscosity. We have studied the case in 
which the plasma was dilute. In this case the coupling 
between the ions is weak. We have found that the vis- 
cosity is constant at low densities and then becomes an 
increasing function of density, contrary to the prediction 
including only binary collisions, where the viscosity re- 
mains constant with density. This fact is due to the cor- 
relations between ions. The dependence of the plasma 
viscosity with temperature, at a fixed density, is the one 
expected for a gas. However, as the temperature de- 
creases changes in the viscosity with density, are larger 
due to stronger ion correlations. 

In the dilute gas limit and for light ion systems, cal- 
culations based on the Chapman-Enskog formalism, in- 
dicate that viscosities with quantum transport cross sec- 
tions are nearly a factor of ten larger than the classical 
results for a fixed screening length of 10 fm. However, 
for larger more realistic screening lengths, the difference 
between quantum and classical calculations is smaller. 

Using the structure factor of the ions we calculate the 
contribution to the shear viscosity due to electron-ion 
scattering. We find that the contribution to the shear 
viscosity from the ions is negligible compared to the for- 
mer one for the conditions of interest in the outer crust of 
neutron stars. Therefore, we do not expect major damp- 
ing of the r-modes from the ions in the crust. 

Our method of calculating the electron contribution 
via the structure factor of the ions also allows us to 
use multi-component plasmas, and the contribution of 
electron-impurity scattering will be automatically in- 
cluded. On the other hand, we have used ion autocor- 
relation functions of the pressure tensor, and the Kubo 
formalism to calculate the viscosity of the ions. This 
method could be extended to calculate other mechanical 
properties of the ions, like the shear modulus, of great 
importance in the understanding of starquakes. 
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Appendix. Calculation of the phase shifts 

Here, we employ an efficient algorithm proposed by 
Klozenberg [5T] to calculate the phase shifts. The radial 
part of the Shrodinger equation is 



u"{r) 



fc^-^-|y(r))..W = 0, (28) 



where /i is the reduced mass, r is the separation distance, 
I is the angular momentum quantum number, k is the 
wave number such that energy E = ^^^ and V{r) is 
a spherically symmetric potential. Two new functions 
ai{r) and Si{r) are defined through the relation 



ui{r) = ai{r) {ji{kr) + Si{r)ni{kr)) 



(29) 



such that 



J/(fcO^^ + Mkr)^^ {Si{r)ai{r)) - 0, (30) 

where ji and ni are the Ricatti-Bessel functions, the so- 
lutions of the free radial equation (28 1 when V{r) = 0. 



With the function Si, Eq. (28) becomes 

^ + ^(.KfcO + ^Kr)n.(fcOf = 0. (31) 
Let us introduce the dimensionless variable t = kr and 



write Eq. (31) as a phase equation 



"^ + Vk{t) im + Si{t)ni{t)f = 0, (32) 

at 



to avoid large values of Si and its infinities. To treat the 
singularities associated with the Ricatti-Bessel functions 
ni(t) when t —^ 0, the normalized analogues 



j,{t)^-10l,niit)^ni{t)t^ 
_ Slit) 



(35) 



21+1 



Slit) ^^, Slit) ^Si{t)t 
are used with the corresponding equations 



dSiit) 21 + 1 
dt ^ 



s 



'lit) + tvkit) (jiit) + Siit)niit)'^ =0, 



(36) 



dSiit) 21 + l~^^^^ _^^^^^^ [Siit)jiit) + hiit))' - 0. 



dt 



t 



(37) 
To find the initial conditions, let us expand the solution 



of Eq. ( 36 1 around t = as 

Siit) = ca + cit + C2t'^ + Oit^), 



(38) 



and assume that the expansion for w(t) is known and has 
the form 



Vkit) ="Y-+ao + ait + Oit^). 



This leads to the coefficients 
Co = 0, 



C2 



Cl 
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(41) 
whence 



Siit<^ 



) ~ Cit + C2t^, 



(42) 



which is the initial condition for Eq. ( 36 I . During numer- 
ical integration, we switch to Eq. (37) when [Sj(i)| > 1, 



and if |5;(t)| > 1 we return back to Eq. (36). This 



2tiV{t/k) 



procedure is carried out until t reaches unity, when we 



where v^it) = ^^p^. By the construction in Eq. (|29|, ^^-^^^ ^^ ^^^ ^^^ ^^ Eq^^ (^ and dMl). After the solu 



the tangent of the phase shift is the limiting value of 
Slit) 



lim Slit) ^t&niSiik)), 

t — '■OO 



(33) 



where the existence of the limit is assured by the fact that 
the potential Vir) falls sufficiently rapidly as r -^ oo. 
When the phase shift is close to | plus multiples of vr, we 
consider the function Siit) — 1/5; (t) and its differential 
equation 



dSiit) 
dt 



Vkit) [Siit)jiit) + niit)y ^ (34) 



tion reaches its limiting value of tan(^/c)) in Eq. (33) to 
satisfactory precision, the evaluation is stopped. 



JWKB phase shifts 

The JWKB approximation allows fast and simple cal- 
culations of phase shifts for large k— and moderately 
large Z-values. Here we briefiy describe this method fol- 
lowing Mott and Massey [32], and, Cohen [53]. Let us 
write the wave equation as 



u'i'ir) + Fir)uiir) = 0, 



(43) 
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A = Z = 1 
a = 10 fm 



and by Jeffrey's method, the solution is 

_ 1 
ui{r) « F]^ " sin 

where now i^i(ro) — 0. 



F^^dr 



(49) 



FIG. 14: (Color online.) Exact and JWKB phase shifts for 
the Yukawa potential in Eq. Q for a system with Z = A = 1 
and o = 10 fm. 



Fir) 



2/i 



(E-Vir)) 
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(44) 



If the particle's energy and angular nionientuni are suf- 
ficiently large, such that the potential energy varies only 
slightly over a few wavelengths, then in the limit h ^ 
we may suppose that F is large. Then the solutions of 
Eq. (43 1 are approximately 



ui 



F 4 exp 



±i / F'sdr 



(45) 



Jeffreys [53] has shown that the solution which decreases 
exponentially inside the classically accessible region (r < 
To) is 



ui[r) fa F 1- sin 



F^-dr 



(46) 



where rp is the classical turning point given by the pos- 
itive root of F(ro) =0. To find the solution which van- 
ishes as r ^ 0, Langer |25] employed the substitutions 



p = log{r) and gi 



_ "i('') ,v 



^/^ 



in Eq. ( 43 1 and obtained 
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For large r, the asymptotic form of Eq. (49 1 is 



+ / {F^ - k)dr + k{r - vq) 



whence 
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■n It: r°° 1 

<5/=^ + y-fcro + y (F/-fc)dr. (51) 



Finally, following Cohen [23], we rewrite the above equa- 
tion in the form used for numerical evaluations: 
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dr 



(52) 



h = ^^^*T^ is the classical impact parameter. 



where ^ — ' 

For numerical computations, we choose two sets of pa- 
rameters for the Yukawa potential in Eq. (B : Z — A — 1 
and a = 10 fm, and, Z = 29 A, A = 88 and a = 26.1 fm. 
An exact calculation of the phase shifts is performed for 
Z = — 5 in the range x = 0~ 1000. For larger values of I, 
up to I = 100 for A = 1 and up to / ^ 500 for A = 88, the 
phase shifts are computed using the the JWKB approxi- 
mation. As shown in Fig. [14] the JWKB results are very 
close to the exact values, even for / = 1. The agreement 
is better for the system with ^ = 88 (not shown here), 
which is essentially classical. 
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